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Abstract 

We propose three new phases of H2O under ultrahigh pressure. Our structural search was 
performed using an adaptive genetic algorithm which allows an extensive exploration of crystal 
structure. The new sequence of pressure-induced transitions beyond ice X at K should be ice 
X— > Pbcra — > Pbca — > Pmc2± — > P2\ — > P2\/c phases. Across the Pmd2\-P2\ transition, the 
coordination number of oxygen increases from 4 to 5 with a significant increase of density. All 
stable crystalline phases have nonmetallic band structures up to 7 TPa. 
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1 



H2O ice is one of the most abundant planet forming materials and its phase diagram is 
of fundamental scientific interest. It is crucial for modeling the interiors of icy solar giants 
(Uranus and Netpune), icy satellites, and also new ocean planets being discovered right now. 
Up to now, sixteen crystalline phases have been identified experimentally P-S]. Ice X is the 
highest-pressure phase among those identified experimentally. In this ionic phase, which 
exists above ~ 70 GPa, oxygen atoms form a bcc lattice and hydrogen atoms are located at 
the midpoint between two nearest-neighbor oxygen atoms. However, in Uranus and Neptune, 
pressure at their core-envelope boundary is estimated to be as high as 0.8 TPa [S], Very 
recently, Neptune-sized icy exoplanets have been discovered [5J [7J- Transitions beyond ice 
X can occur under these high pressure conditions and they are essential for modeling the 
interiors of these planets. Previously, Benoit et al. [8j and Caracas |9] predicted a phase 
with Pbcm symmetry as the next high-pressure phase after ice X (up to 0.3 TPa). More 
recently, Militzer and Wilson proposed transitions from the Pbcm to Pbca and Cmcm phases 
at 0.76 and 1.55 TPa, respectively [10J. Interestingly, they predicted the Cmcm phase to be 
metallic, an important property for understanding the origin of magnetic fields in the giants. 
While these studies surely predict the existence of new phases, their searches using MD and 
phonon instability calculations explored limited regions of phase space producing structure 
relatively close to that of ice X. In this Letter, we report new structures of solid H2O in the 
TPa pressure regime. They were found using a novel, general, and efficient structural search 
algorithm, the adaptive genetic algorithm. 

The search for the lowest-enthalpy structures of ultrahigh-pressure ice is based on the 
Deaven-Ho genetic algorithm (GA) scheme jllj . This method combined with first-principles 
calculation has been proven to work very efficiently [T2Tll5| . However, first-principles cal- 
culations are computationally very demanding for the GA scheme when there are a large 
number of atoms. On the other hand, GA searches using empirical model potentials are fast 
but suffer from inaccuracies which can lead the search to wrong structures. However, most 
of the child structures generated in GA are actually not favorable except in the simplest 
problems; many false structures have to be tried before hitting on the correct structure. 
Most of the computer time is spent in relaxing a large number of false candidate structures. 
To remedy this, we introduced a new structure search technique, the adaptive GA. In the 
new scheme, we employ auxiliary model potentials to estimate the energy ordering of dif- 
ferent competing geometries in a preliminary stage. Parameters of the auxiliary potentials 
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are adaptively adjusted to reproduce first-principles results during the course of the GA 
search. This adaptive approach also allows the system to hop from one basin to another in 
the energy landscape leading to efficient sampling of configuration space. 

In our study, we found that the packing geometry, volumes(pressures) and energy ordering 
of ice crystal phases at high pressures (>0.5 TPa) can be described relatively well by simple 
Lennard- Jones (LJ) potentials: 



In the H2O system, a total of six parameters, representing the O-O, 0-H and H-H inter- 
actions, are adaptively adjusted to explore structural phase space during the GA search. 
We initiate the search with random structures and carry out DFT calculations to get their 
energies, atomic forces and stress tensor. The LJ potential parameters are fitted to these 
structure by force-matching method[16j. Based on this auxiliary potential we perform full 
GA search to yield new structural candidates. DFT calculations are carried out on the new 
GA pool and potential parameters are adaptively adjusted again. These procedures can 
be done iteratively until the potential parameters converge. The final LJ potential pool 
population will then be examined with full DFT relaxation. This adaptive scheme allows a 
very efficient sampling of crystal structures, we can easily search relatively large systems up 
to 12 H2O formula units. 

Details of the GA search can be found in several previous papers |T2HT5] . For auxiliary 
potential optimization we use LAMMPS code [17] and conjugate gradient method. First 
principles calculations are performed using density-functional theory (DFT) within PBE- 
type generalized-gradient approximation (GGA) [18j. Vanderbilt-type pseudopotentials |19| 
were generated using the following valence electron configurations: Is 1 and 2s 2 2p 4 with 
cutoff radii of 0.5 and 1.4 a.u for hydrogen and oxygen, respectively. Candidate structures 
obtained with the adaptive GA were then refined using a harder oxygen pseudopotential with 
valence electron configuration of 2s 2 2p 4 3d° and cutoff radius of 1.0 a.u. The cutoff energy 
for plane-wave expansions are 40 Ry and 120 Ry for the softer and the harder pseudopo- 
tentials respectively. Brillouin-zone integration was performed using the Monkhorst-Pack 
sampling scheme [20] over k-point meshes of spacing 2tt x 0.05A -1 . Structural relaxations 
were performed using variable-cell-shape molecular dynamics [2T1 [22] . To test for structural 
stability, phonon and vibrational density of states (VDOS) calculations were carried out 
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for candidate structures using density-functional-perturbation theory [23l 124"] . Zero-point 
motion (ZPM) effects are taken into account within the quasiharmonic approximation [25J. 
All first-principles calculations were performed using the Quantum-ESPRESSO [26], which 
has been interfaced with the GA scheme in a fully parallel manner. 

Fig. [T] shows average values of the pressures and enthalpies calculated by first principles 
for structures in the adaptive GA pool. Target pressure of this adaptive GA search is 2 TPa. 
Our results show a fast convergence of the adaptive GA. After 10 iterations, LJ-potential 
pressures for structures in the GA pool are almost identical to first-principles DFT results. 
The adaptive GA search successfully predicts three new structures: Pmc2\ at 1 TPa, P2\ 
at 2 TPa, and P2\/c phase at 3 TPa (Fig. [2j. They consist of 4, 4, and 8 formula units, 
respectively. Structures with 3, 5, 6, 7, 9, 10, and 11 formula units were not energetically 
competitive. Pbcm and Pbca phases proposed previously were also obtained as metastable 
phases at these pressures. Examination of the enthalpy differences between the different 
competing phases (Fig. [3j shows that the sequence of pressure-induced phase transitions 
beyond ice X is X— >■ Pbcm — > Pbca — > Pmc2\ — > P2i — > P2\jc. Static transition pressures 
are 0.28, 0.75, 0.89, 1.28, and 2.68 TPa, respectively. Zero-point motion strongly affects 
transition pressures. For X-Pbcm, Pbcm-Pbca, and Pbca-Pmc2\ transitions, ZPM increases 
transition pressures to 0.29, 0.77, and 0.92 TPa. On the other hand, for Pmc2\-P2\ and 
P2 1 -P2 1 /c transitions, ZPM greatly decreases transition pressures to 1.14 and 1.96 TPa. 

The structures of these competing phases are closely related. Under compression, ice X 
transforms to the Pbcm phase and then to the Pbca phase by means of soft phonon related 
deformations |9], HO]. In these three phases, all hydrogen atoms are located at the midpoint 
between two neighboring oxygen atoms. However, during the transition to the Pmc2\ phase, 
two among the four hydrogen atoms bonded to the 01 oxygen jump to the octahedral 
interstitial sites next to two second-nearest-neighbor oxygen atoms. The arrangements of 
hydrogen atoms are Pbca-like around 01 and Cmcm-like around 02. Therefore, this phase 
is structurally intermediate between the Pbca and Cmcm phases. Its becomes metastable 
with respect to the Cmcm phase at ~2.5 TPa. In ice X, Pbcm, Pbca, and Pmc2\ phases, 
OH4 tetrahedra form the basic structural unit, with connectivity varying in each phase. 
Across the Pbca to Pmc2\ transition, two interpenetrating tetrahedral networks transform 
into a single network. Locally, hydrogen atoms around 02 try to keep two interpenetrating 
networks, while those around 01 atoms connect two networks into a single one. Tetrahedra 
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around the 01 atoms are severely distorted compared to those around the 02 atoms. 

The Pmc2i phase is dynamically stable up to 1.3 TPa. At ~1.5 TPa, a zone-center soft 
mode appears. This soft mode induces a monoclinic distortion giving rise to the P2\ phase. 
The P2\ space group is a subgroup of Pmc2\. In the P2\ phase, the structural unit is no 
longer OH4. The coordination number of oxygen atoms increases from 4 to 5. In X, Pbcm, 
Pbca, and Pmc2\ phases, phonon dispersions are divided into two groups corresponding 
to nearly-rigid-body and internal 0-H stretching motions of 0H 4 tetrahedra, while this 
distinction is blurred in the P2\ phase. As a result of higher coordination number, increase 
in H-0 bond-length, and a 2.0% volume reduction across the Pmc2 1 -P2 1 transition. Both 
P2\ and P2\jc phases are dynamically stable at least up to 7 TPa, with no soft modes 
developing under compression. 

In contrast with the metallic Cmcm phase predicted by |10j . all three new phases (Pmc2i, 
P2i, P2i/c) have substantial DFT band gaps. In the P2i and P2\jc phases, the band 
gap decreases slowly under compression (Fig. [4]), closing at ~7 TPa in the P2\jc phase. 
Although all crystalline phases are insulating, H 2 could be a good conductor at relevant 
conditions for two reasons: 1) protons highly mobile in the oxygen sublattice producing 
superionic phases at ultrahigh pressures and temperatures typical of the interiors of icy 
solar giants and exoplanets; 2) depending on the temperature, thermally excited carriers 
also contribute to increase the conductivity [27H34| . So far, only the bcc oxygen sub-lattice 
has been considered in the investigation of superionic phases. Our new crystal structures 
provide starting points for further investigation of conducting superionic states of these 
phases at ultrahigh pressures and temperatures. They indicate that the oxygen sublattice 
should prefer to have hexagonal-derived structures beyond ~0.4 TPa. This possibly implies 
a bcc-hcp transition in the superionic phase. 
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Figure 1: Convergence of average DFT pressure and enthalpy of structures in the LJ-potential 
GA pool as a function of adaptive potential iterations. The DFT pressure converges to the target 
pressure of 2 TPa after 10 iterations. The unit cell contains 8 H2O units. 



Pbcm Pbca Pmc2 1 P2 1 P2 1 /c 




Figure 2: Crystal structures of ultrahigh-pressure phases of ice. Blue and red large spheres denote 
oxygen atoms. White small spheres denote hydrogen atoms. 
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Pmc2i-type H 2 at 1 TPa 
(a,b,c) (3.087A, 1.890A, 3.296A) 



Hi 


4c 


(0.25235, 0.39936, 0.34521) 


H 2 


2a 


(0, 0.72441, 0.58937) 


TT 

H 3 


2ft 


(0.5, 0.02194, 0.05554) 


r\ 

Ui 


2a 


(0, 0.79553, -0.00013) 




2ft 


(0.5, 0.71235, 0.29490) 






P2i-type H 2 at 2 TPa 


(a,b, 


c,/3) (1.711A, 3.066A, 2.825A, 99.83°) 


Hi 


2a 


(0.03146, -0.004190, 0.97579) 


H 2 


2a 


(0.17734, 0.60082, 0.33256) 


H 3 


2a 


(0.25493, 0.38266, 0.73104) 


IT 


2a 


(0.56268, 0.74399, 0.68952) 


Oi 


2a 


(0.82524, 0.52141, 0.52159) 




2a 


(0.34462, 0.75497, 0.01421) 






P2i/c-type H 2 at 3 TPa 


(a, ft, 


c,p) (2.921A, 2.890A, 3.338A, 117.86°) 


Hi 


4e 


(-0.49599, 0.19043, 0.44690) 


H 2 


4e 


(-0.14264, 0.13088, -0.07153) 


H 3 


4e 


(-0.24986, -0.50964, -0.27537) 


H 4 


4e 


(0.21275, 0.37363, -0.03277) 


Oi 


4e 


(0.06825, -0.36013, -0.15960) 


o 2 


4e 


(-0.42335, -0.37490, 0.33745) 



Table I: Structural parameters of Pmc2\-, P2\-, and P2i/c-type H 2 0. 
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Figure 3: Enthalpies of ultrahigh-pressure phases of ice. Dashed vertical lines denote static tran- 
sition pressures. 
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Figure 4: Electronic band gaps of ultrahigh-pressure phases of ice. Dashed vertical lines denote 
static transition pressures. 
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